##### 加载包、自定义函数
# 数据输入输出
library(readxl)    # 读取速度优于 read.xlsx()，且不依赖 Java
library(xlsx)      # 输出数据时要用到 write.xlsx()
library(openxlsx)  # 效率读、写 xlsx，以及定制 xlsx
library(officer)
library(flextable)
library(gtsummary)
# 数据处理
library(reshape2)
# 分析
library(geepack)
library(mice)

# 加载自定义函数
source("./gee_run.r")
#####


# 5到10行不用run!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
t1 = read_excel("./data918.xlsx", sheet = 3)
t2 = read_excel("./data918.xlsx", sheet = 4)
extract_id = unique(t1$id)
data = merge(t1,t2,by="id",all=TRUE)
data_final = data[data$id %in% extract_id, ]
write.xlsx(data_final, './merge.xlsx', overwrite = T)
# 从这下面开始读
t3 = read_excel("./merge.xlsx", sheet = 1)
data_test = t3
# data_test = data_test[data_test$score_class12 == 1,]  # 只取了12月不合格的对象


# 1.1 删除 6 月的观测、且 8 月只有一个观测的个体
data_cast = dcast(data_test, id~month)
data_cast = transform(data_cast, count=`8`+`12`+`18`+`24`)
extract_id2 = data_cast[data_cast$count>1, ]$id
data_tmp = data_test[(data_test$id %in% extract_id2) & (data_test$month != 6), ]
# 1.2 提取数据，确定自变量、因变量、协变量
data_anal = data_tmp[, c(1:2, 9:13, 14:25, 27, 29:39,63,65)]
varsx = colnames(data_test)[c(14:25, 65)]  # 自变量
varsy = colnames(data_test)[c(9:13)]       # 因变量
varsc = colnames(data_test)[c(29:39)]      # 协变量
# 1.3 确保变量类型一致：因变量为数值，自变量、协变量为因子
data_anal$id = as.character(data_anal$id)
data_anal[varsy] = lapply(data_anal[varsy], as.numeric)
data_anal[varsx] = lapply(data_anal[varsx], as.factor)
data_anal[varsc] = lapply(data_anal[varsc], as.factor)


# 2. 填补缺失值并跑 GEE
# 2.1 mice
data_tmp = unique(data_anal[, c("id", varsc)])
miss = md.pattern(data_tmp)                # 查看缺失数据的模式
imp = mice(data_tmp, m = 5, seed = 1)      # 填补
mice_data = complete(imp, action = 1)      # 懒惰的人选择第一个
mice_data = merge(data_anal[, c("id", "month", varsy, varsx)], mice_data)   # 合并
write.xlsx(mice_data, './AVERstratified_mice_data.xlsx', overwrite = T)

## 2.2 预处理
mice_data = read.xlsx('./AVERstratified_mice_data.xlsx', 1)
mice_data = mice_data[complete.cases(mice_data[, varsx]), ]   # 删除自变量为空的行
mice_data$score_qualified = ifelse(mice_data$score_qualified == 1, 2, 1)  # 置换  score_qualified（1 -> 2, 2 -> 1）
mice_data$id = as.character(mice_data$id)
mice_data[varsy] = lapply(mice_data[varsy], as.numeric)
mice_data[varsx] = lapply(mice_data[varsx], as.factor)
mice_data[varsc] = lapply(mice_data[varsc], as.factor)

# 2.3 gee
fit_t1 = gee_run(
  mice_data,
  varsx = varsx,
  varsy = varsy,
  varsc = varsc[-5],
  file = "./AVERstratified_10cv_del_birthweight_c.xlsx"
)
fit_t2 = gee_run(
  mice_data,
  varsx = varsx,
  varsy = varsy,
  varsc = varsc[-2],
  file = "../output/1009_mice_10cv_del_week_c.xlsx"
)



# 加载自定义函数
# give a fit_table, return a summary_plot_table
summary_out = function (fit_table) {
  # summary_out
  ret = list()
  for (i in 1:length(fit_table[[1]])) {
    tb1 = list()
    for (j in 1:length(fit_table)) {
      tb1[[j]] = tbl_regression(fit_table[[j]][[i]], exponentiate = T, include = c(1))
      print(i)
    }
    ret[[i]] = tbl_merge(tb1, tab_spanner = varsy)
  }
  return (ret)
}

s1 = summary_out(fit_t1)
sect_properties = prop_section(
  page_size = page_size(orient = "landscape", width = 8.3, height = 11.7),
  type = "continuous",
  page_margins = page_mar()
)
save_as_docx(
  values = lapply(s1, as_flex_table),
  path = "./output/921_summary_delete.docx", 
  pr_section = sect_properties
)
